#==========================================================================
# Mauerer, I. & Puy, M. S. (2025). 
# "Ingroup and Outgroup Effects on Party Placement Perceptions" 
# European Political Science Review
#==========================================================================

#---------------------------------------------------------------------------------------
# This file contains the commands to run the models and create the tables and figures
# R version 4.5.1 
# Platform: x86_64-w64-mingw32/x64 (64-bit)
#----------------------------------------------------------------------------------------


#==================
# Packages
#==================
library(xtable)
library(ggplot2)
require(gridExtra)
library(data.table)

#==================
# Working Directory
#==================
setwd("...")


#==================
# Load data
#==================
load("Data.RDS")

#===========================
# Define reference groups
#===========================
basque20$natid  <- relevel(basque20$natid , ref = "2. neutral")
basque20$status <- relevel(basque20$status, ref = "2. middle")
basque20$ideo5  <- relevel(basque20$ideo5 , ref = "Nationalist")

basque16$natid  <- relevel(basque16$natid , ref = "2. neutral")
basque16$status <- relevel(basque16$status, ref = "2. middle")
basque16$ideo5  <- relevel(basque16$ideo5 , ref = "Nationalist")

basque12$natid  <- relevel(basque12$natid , ref = "2. neutral")
basque12$status <- relevel(basque12$status, ref = "2. middle")
basque12$ideo5  <- relevel(basque12$ideo5 , ref = "Nationalist")

basque09$natid  <- relevel(basque09$natid , ref = "2. neutral")
basque09$status <- relevel(basque09$status, ref = "2. middle")
basque09$ideo5  <- relevel(basque09$ideo5 , ref = "Nationalist")


#====================================================
# Table 1: Mean Left-Right Self-Placements by Groups
#====================================================

#---------------------
# Left-Right [lr_ego]
#---------------------
means_lrego <- c(round(mean(basque20$lr_ego), digits=3), 
                 round(mean(basque16$lr_ego), digits=3),
                 round(mean(basque12$lr_ego), digits=3),
                 round(mean(basque09$lr_ego), digits=3))
#
sd_lrego    <- c(round(sd(basque20$lr_ego), digits=3), 
                 round(sd(basque16$lr_ego), digits=3),
                 round(sd(basque12$lr_ego), digits=3),
                 round(sd(basque09$lr_ego), digits=3))

#--------------------------
# National Groups [natid]
#--------------------------
means_lrego_natid20   <- tapply(basque20$lr_ego, basque20$natid,mean)
means_lrego_natid16   <- tapply(basque16$lr_ego, basque16$natid,mean)
means_lrego_natid12   <- tapply(basque12$lr_ego, basque12$natid,mean)
means_lrego_natid09   <- tapply(basque09$lr_ego, basque09$natid,mean)
#
sd_lrego_natid20      <- tapply(basque20$lr_ego, basque20$natid,sd)
sd_lrego_natid16      <- tapply(basque16$lr_ego, basque16$natid,sd)
sd_lrego_natid12      <- tapply(basque12$lr_ego, basque12$natid,sd)
sd_lrego_natid09      <- tapply(basque09$lr_ego, basque09$natid,sd)

# put years together
means_lrego_natid <- cbind(round(means_lrego_natid20, digits=3), 
                           round(means_lrego_natid16, digits=3), 
                           round(means_lrego_natid12, digits=3),
                           round(means_lrego_natid09, digits=3))
#
sd_lrego_natid    <- cbind(round(sd_lrego_natid20, digits=3), 
                           round(sd_lrego_natid16, digits=3), 
                           round(sd_lrego_natid12, digits=3),
                           round(sd_lrego_natid09, digits=3))

#-----------------------------
# Ideological Groups [ideo5]
#-----------------------------
means_lrego_ideo520   <- tapply(basque20$lr_ego, basque20$ideo5,mean)
means_lrego_ideo516   <- tapply(basque16$lr_ego, basque16$ideo5,mean)
means_lrego_ideo512   <- tapply(basque12$lr_ego, basque12$ideo5,mean)
means_lrego_ideo509   <- tapply(basque09$lr_ego, basque09$ideo5,mean)
#
sd_lrego_ideo520      <- tapply(basque20$lr_ego, basque20$ideo5,sd)
sd_lrego_ideo516      <- tapply(basque16$lr_ego, basque16$ideo5,sd)
sd_lrego_ideo512      <- tapply(basque12$lr_ego, basque12$ideo5,sd)
sd_lrego_ideo509      <- tapply(basque09$lr_ego, basque09$ideo5,sd)

# put years together
means_lrego_ideo5 <- cbind(round(means_lrego_ideo520, digits=3),  
                           round(means_lrego_ideo516, digits=3),  
                           round(means_lrego_ideo512, digits=3),  
                           round(means_lrego_ideo509, digits=3))
#
sd_lrego_ideo5    <- cbind(round(sd_lrego_ideo520, digits=3),  
                           round(sd_lrego_ideo516, digits=3),  
                           round(sd_lrego_ideo512, digits=3),  
                           round(sd_lrego_ideo509, digits=3))

#---------
# gender
#---------
means_lrego_gender20   <- tapply(basque20$lr_ego, basque20$gender,mean)
means_lrego_gender16   <- tapply(basque16$lr_ego, basque16$gender,mean)
means_lrego_gender12   <- tapply(basque12$lr_ego, basque12$gender,mean)
means_lrego_gender09   <- tapply(basque09$lr_ego, basque09$gender,mean)
#
sd_lrego_gender20      <- tapply(basque20$lr_ego, basque20$gender,sd)
sd_lrego_gender16      <- tapply(basque16$lr_ego, basque16$gender,sd)
sd_lrego_gender12      <- tapply(basque12$lr_ego, basque12$gender,sd)
sd_lrego_gender09      <- tapply(basque09$lr_ego, basque09$gender,sd)

# put years together
means_lrego_gender <- cbind(round(means_lrego_gender20, digits=3),  
                            round(means_lrego_gender16, digits=3),  
                            round(means_lrego_gender12, digits=3),  
                            round(means_lrego_gender09, digits=3))
#
sd_lrego_gender    <- cbind(round(sd_lrego_gender20, digits=3),  
                            round(sd_lrego_gender16, digits=3),  
                            round(sd_lrego_gender12, digits=3),  
                            round(sd_lrego_gender09, digits=3))

#-----------------------------------------
# Catholic [cath] (not contained in 2009) 
#-----------------------------------------
means_lrego_cath20   <- tapply(basque20$lr_ego, basque20$cath,mean)
means_lrego_cath16   <- tapply(basque16$lr_ego, basque16$cath,mean)
means_lrego_cath12   <- tapply(basque12$lr_ego, basque12$cath,mean)
#
sd_lrego_cath20      <- tapply(basque20$lr_ego, basque20$cath,sd)
sd_lrego_cath16      <- tapply(basque16$lr_ego, basque16$cath,sd)
sd_lrego_cath12      <- tapply(basque12$lr_ego, basque12$cath,sd)

# put years together
means_lrego_cath  <- cbind(round(means_lrego_cath20, digits=3),  
                           round(means_lrego_cath16, digits=3),  
                           round(means_lrego_cath12, digits=3),  NA)
# 
sd_lrego_cath     <- cbind(round(sd_lrego_cath20, digits=3),  
                           round(sd_lrego_cath16, digits=3),  
                           round(sd_lrego_cath12, digits=3),  NA)


#------------------------
# Social class [status]
#------------------------
means_lrego_status20   <- tapply(basque20$lr_ego, basque20$status,mean)
means_lrego_status16   <- tapply(basque16$lr_ego, basque16$status,mean)
means_lrego_status12   <- tapply(basque12$lr_ego, basque12$status,mean)
means_lrego_status09   <- tapply(basque09$lr_ego, basque09$status,mean)
#
sd_lrego_status20      <- tapply(basque20$lr_ego, basque20$status,sd)
sd_lrego_status16      <- tapply(basque16$lr_ego, basque16$status,sd)
sd_lrego_status12      <- tapply(basque12$lr_ego, basque12$status,sd)
sd_lrego_status09      <- tapply(basque09$lr_ego, basque09$status,sd)

# put years together
means_lrego_status <- cbind(round(means_lrego_status20, digits=3),  
                            round(means_lrego_status16, digits=3),  
                            round(means_lrego_status12, digits=3),  
                            round(means_lrego_status09, digits=3))
#
sd_lrego_status    <- cbind(round(sd_lrego_status20, digits=3),  
                            round(sd_lrego_status16, digits=3),  
                            round(sd_lrego_status12, digits=3),  
                            round(sd_lrego_status09, digits=3))

#---------------------
# put groups together
#---------------------
table1 <- rbind(means_lrego,        sd_lrego,
                means_lrego_natid,  sd_lrego_natid,
                means_lrego_ideo5,  sd_lrego_ideo5,
                means_lrego_gender, sd_lrego_gender, 
                means_lrego_cath,   sd_lrego_cath, 
                means_lrego_status, sd_lrego_status)
colnames(table1) <- c(2020, 2016, 2012, 2009)
table1



#=========================================
# Multivariate Linear Model Estimations
#=========================================

#----------------
# 2020 Election
#----------------
m20 <- lm(cbind(lr_podemos,lr_pse, lr_pp, lr_eh, lr_pnv) ~ natid + eusk + lr_ego + ideo5 + eco + gender + aged + status + cath + pop, 
          data=basque20)
m20
summary(m20)

#----------------
# 2016 Election
#----------------
m16 <- lm(cbind(lr_podemos, lr_pse, lr_pp, lr_eh, lr_pnv) ~ natid + eusk + lr_ego + ideo5 + eco + gender + aged + status + cath + pop, 
          data=basque16)
m16
summary(m16)

#----------------
# 2012 Election 
#----------------
m12 <- lm(cbind(lr_iu, lr_pse, lr_pp, lr_eh, lr_pnv) ~ natid + eusk12 + lr_ego + ideo5 + eco + gender + aged + status + cath + pop, 
          data=basque12)
m12
summary(m12)

#----------------
# 2009 Election
#----------------
m09 <- lm(cbind(lr_eb,  lr_pse, lr_pp,lr_aralar, lr_pnv) ~ natid + eusk + lr_ego + ideo5 +  eco + gender + aged + status + pop, 
          data=basque09)
m09
summary(m09)



#====================================================
# Table 2: Persistent Perceptual Group Divergences
#====================================================

# (a) National Groups
table2a <- matrix(NA, 8, ncol = 5)
rownames(table2a)<- c(2020, 2016, 2012, 2009, 2020, 2016, 2012, 2009)
colnames(table2a)<- c("PSE", "PP", "Podemos", "PNV", "EH-Bildu")


# Basque - Neutral Comparison 
table2a[1:4, 1] <- round(rbind(summary(m20)$`Response lr_pse`$coefficients[3,1],
                               summary(m16)$`Response lr_pse`$coefficients[3,1],
                               summary(m12)$`Response lr_pse`$coefficients[3,1],
                               summary(m09)$`Response lr_pse`$coefficients[3,1]), digits=3)

table2a[1:4, 2] <- round(rbind(summary(m20)$`Response lr_pp`$coefficients[3,1],
                               summary(m16)$`Response lr_pp`$coefficients[3,1],
                               summary(m12)$`Response lr_pp`$coefficients[3,1],
                               summary(m09)$`Response lr_pp`$coefficients[3,1]), digits=3)

table2a[1:4, 3] <- round(rbind(summary(m20)$`Response lr_podemos`$coefficients[3,1],
                               summary(m16)$`Response lr_podemos`$coefficients[3,1],
                               NA,
                               summary(m09)$`Response lr_eb`$coefficients[3,1]), digits=3)
 
# Spanish - Neutral Comparison 
table2a[5:8, 2] <- round(rbind(summary(m20)$`Response lr_pp`$coefficients[2,1],
                               summary(m16)$`Response lr_pp`$coefficients[2,1],
                               NA,
                               summary(m09)$`Response lr_pp`$coefficients[2,1]), digits=3)

table2a



# (b) Ideological Groups
table2b <- matrix(NA, 16, ncol = 5)
rownames(table2b)<- c(2020, 2016, 2012, 2009, 2020, 2016, 2012, 2009,
                      2020, 2016, 2012, 2009, 2020, 2016, 2012, 2009)
colnames(table2b)<- c("PSE", "PP", "Podemos", "PNV", "EH-Bildu")


# Conservative - Nationalist Comparison
table2b[1:4, 3] <- round(rbind(summary(m20)$`Response lr_podemos`$coefficients[6,1],
                               summary(m16)$`Response lr_podemos`$coefficients[6,1],
                               summary(m12)$`Response lr_iu`$coefficients[6,1],
                               NA), digits=3)

table2b[1:4, 4] <- round(rbind(NA,
                               summary(m16)$`Response lr_pnv`$coefficients[6,1],
                               summary(m12)$`Response lr_pnv`$coefficients[6,1],
                               summary(m09)$`Response lr_pnv`$coefficients[6,1]), digits=3)

#Socialist - Nationalist Comparison
table2b[5:8, 4] <- round(rbind(summary(m20)$`Response lr_pnv`$coefficients[7,1],
                               summary(m16)$`Response lr_pnv`$coefficients[7,1],
                               summary(m12)$`Response lr_pnv`$coefficients[7,1],
                               summary(m09)$`Response lr_pnv`$coefficients[7,1]), digits=3)

# Liberal - Nationalist Comparison
table2b[9:12, 2] <- round(rbind(summary(m20)$`Response lr_pp`$coefficients[9,1],
                               summary(m16)$`Response lr_pp`$coefficients[9,1],
                               NA,
                               summary(m09)$`Response lr_pp`$coefficients[9,1]), digits=3)

# Feminist/Ecologist - Nationalist Comparison
table2b[13:16, 4] <- round(rbind(NA,
                              summary(m16)$`Response lr_pnv`$coefficients[8,1],
                              summary(m12)$`Response lr_pnv`$coefficients[8,1],
                              summary(m09)$`Response lr_pnv`$coefficients[8,1]), digits=3)

table2b

# (c) Socioeconomic Groups
table2c <- matrix(NA, 3, ncol = 5)
rownames(table2c)<- c(2020, 2016, 2012)
colnames(table2c)<- c("PSE", "PP", "Podemos", "PNV", "EH-Bildu")      
      
table2c[1:3, 1] <- round(rbind(summary(m20)$`Response lr_pse`$coefficients[15,1],
                               summary(m16)$`Response lr_pse`$coefficients[15,1],
                               summary(m12)$`Response lr_pse`$coefficients[15,1]), digits=3)

table2c[1:3, 4] <- round(rbind(summary(m20)$`Response lr_pnv`$coefficients[15,1],
                               summary(m16)$`Response lr_pnv`$coefficients[15,1],
                               NA), digits=3)

table2c[1:3, 5] <- round(rbind(summary(m20)$`Response lr_eh`$coefficients[15,1],
                               NA,
                               summary(m12)$`Response lr_eh`$coefficients[15,1]), digits=3)

table2c


#-------------------------------------------------------------
# Online Appendix 
#-------------------------------------------------------------


#-------------------------------------------------------------
# B Mean Left-Right Party Placement Perceptions
#-------------------------------------------------------------

#--------------------------------------------------------------------------
# Table A1: Mean Left-Right Party Placement Perceptions across all Groups
#--------------------------------------------------------------------------

# 2020: 
frameLR20 =data.frame(podemos=basque20$lr_podemos, 
                      pse=basque20$lr_pse, 
                      pp=basque20$lr_pp, 
                      eh=basque20$lr_eh, 
                      pnv=basque20$lr_pnv)

frameLR20T<-data.table(frameLR20)

lr20_party<- rbind(frameLR20T[,list(mean=mean(podemos))],
                   frameLR20T[,list(mean=mean(pse))],
                   frameLR20T[,list(mean=mean(pp))],
                   frameLR20T[,list(mean=mean(eh))],
                   frameLR20T[,list(mean=mean(pnv))])

lr20_party <- round(lr20_party, digits=3)

# 2016: 
frameLR16 =data.frame(podemos=basque16$lr_podemos, 
                      pse=basque16$lr_pse, 
                      pp=basque16$lr_pp, 
                      eh=basque16$lr_eh, 
                      pnv=basque16$lr_pnv)

frameLR16T<-data.table(frameLR16)

lr16_party<- rbind(frameLR16T[,list(mean=mean(podemos))],
                   frameLR16T[,list(mean=mean(pse))],
                   frameLR16T[,list(mean=mean(pp))],
                   frameLR16T[,list(mean=mean(eh))],
                   frameLR16T[,list(mean=mean(pnv))])

lr16_party <- round(lr16_party, digits=3)


# 2012: 
frameLR12 =data.frame(iu=basque12$lr_iu, 
                      pse=basque12$lr_pse, 
                      pp=basque12$lr_pp, 
                      eh=basque12$lr_eh, 
                      pnv=basque12$lr_pnv)

frameLR12T<-data.table(frameLR12)

lr12_party<- rbind(frameLR12T[,list(mean=mean(iu))],
                   frameLR12T[,list(mean=mean(pse))],
                   frameLR12T[,list(mean=mean(pp))],
                   frameLR12T[,list(mean=mean(eh))],
                   frameLR12T[,list(mean=mean(pnv))])

lr12_party <- round(lr12_party, digits=3)


# 2009:
frameLR09 =data.frame(eb=basque09$lr_eb, 
                      pse=basque09$lr_pse, 
                      pp=basque09$lr_pp, 
                      aralar=basque09$lr_aralar, 
                      pnv=basque09$lr_pnv)

frameLR09T<-data.table(frameLR09)

lr09_party<- rbind(frameLR09T[,list(mean=mean(eb))],
                   frameLR09T[,list(mean=mean(pse))],
                   frameLR09T[,list(mean=mean(pp))],
                   frameLR09T[,list(mean=mean(aralar))],
                   frameLR09T[,list(mean=mean(pnv))])

lr09_party <- round(lr09_party, digits=3)


tableA1<- cbind(lr20_party, lr16_party, lr12_party, lr09_party)

colnames(tableA1)<-c("2020", "2016", "2012", "2009")
rownames(tableA1)<-c("Podemos", "PSE", "PP", "EH-Bildu", "PNV") 

tableA1



#------------------------------------------------------------------------------------
# Table A2: Mean Left-Right Party Placement Perceptions by Groups: Podemos
#------------------------------------------------------------------------------------

#----------
# natid
#----------
means_lr_podemos_natid20   <- tapply(basque20$lr_podemos, basque20$natid,mean)
means_lr_podemos_natid16   <- tapply(basque16$lr_podemos, basque16$natid,mean)
means_lr_podemos_natid12   <- tapply(basque12$lr_iu, basque12$natid,mean)
means_lr_podemos_natid09   <- tapply(basque09$lr_eb, basque09$natid,mean)

# put years together
means_lr_podemos_natid <- cbind(round(means_lr_podemos_natid20, digits=3), 
                                round(means_lr_podemos_natid16, digits=3), 
                                round(means_lr_podemos_natid12, digits=3),
                                round(means_lr_podemos_natid09, digits=3))
means_lr_podemos_natid

#---------
# ideo5
#---------
means_lr_podemos_ideo520   <- tapply(basque20$lr_podemos, basque20$ideo5,mean)
means_lr_podemos_ideo516   <- tapply(basque16$lr_podemos, basque16$ideo5,mean)
means_lr_podemos_ideo512   <- tapply(basque12$lr_iu, basque12$ideo5,mean)
means_lr_podemos_ideo509   <- tapply(basque09$lr_eb, basque09$ideo5,mean)

# put years together
means_lr_podemos_ideo5 <- cbind(round(means_lr_podemos_ideo520, digits=3),  
                                round(means_lr_podemos_ideo516, digits=3),  
                                round(means_lr_podemos_ideo512, digits=3),  
                                round(means_lr_podemos_ideo509, digits=3))
means_lr_podemos_ideo5

#---------
# gender
#---------
means_lr_podemos_gender20   <- tapply(basque20$lr_podemos, basque20$gender,mean)
means_lr_podemos_gender16   <- tapply(basque16$lr_podemos, basque16$gender,mean)
means_lr_podemos_gender12   <- tapply(basque12$lr_iu, basque12$gender,mean)
means_lr_podemos_gender09   <- tapply(basque09$lr_eb, basque09$gender,mean)


# put years together
means_lr_podemos_gender <- cbind(round(means_lr_podemos_gender20, digits=3),  
                                 round(means_lr_podemos_gender16, digits=3),  
                                 round(means_lr_podemos_gender12, digits=3),  
                                 round(means_lr_podemos_gender09, digits=3))
means_lr_podemos_gender


#---------
# cath (not contained in 2009)
#---------
means_lr_podemos_cath20   <- tapply(basque20$lr_podemos, basque20$cath,mean)
means_lr_podemos_cath16   <- tapply(basque16$lr_podemos, basque16$cath,mean)
means_lr_podemos_cath12   <- tapply(basque12$lr_iu, basque12$cath,mean)


# put years together
means_lr_podemos_cath <- cbind(round(means_lr_podemos_cath20, digits=3),  
                               round(means_lr_podemos_cath16, digits=3),  
                               round(means_lr_podemos_cath12, digits=3), NA)
means_lr_podemos_cath



#---------
# status
#---------
means_lr_podemos_status20   <- tapply(basque20$lr_podemos, basque20$status,mean)
means_lr_podemos_status16   <- tapply(basque16$lr_podemos, basque16$status,mean)
means_lr_podemos_status12   <- tapply(basque12$lr_iu, basque12$status,mean)
means_lr_podemos_status09   <- tapply(basque09$lr_eb, basque09$status,mean)


# put years together
means_lr_podemos_status <- cbind(round(means_lr_podemos_status20, digits=3),  
                                 round(means_lr_podemos_status16, digits=3),  
                                 round(means_lr_podemos_status12, digits=3),  
                                 round(means_lr_podemos_status09, digits=3))
means_lr_podemos_status


#---------------------
# put variables together
#---------------------
tableA2 <- rbind(means_lr_podemos_natid,means_lr_podemos_ideo5, 
                           means_lr_podemos_gender, means_lr_podemos_cath, means_lr_podemos_status)
colnames(tableA2) <- c(2020, 2016, 2012, 2009)

tableA2 


#------------------------------------------------------------------------------------
# Table A3: Mean Left-Right Party Placement Perceptions by Groups: PSE
#------------------------------------------------------------------------------------

#----------
# natid
#----------
means_lr_pse_natid20   <- tapply(basque20$lr_pse, basque20$natid,mean)
means_lr_pse_natid16   <- tapply(basque16$lr_pse, basque16$natid,mean)
means_lr_pse_natid12   <- tapply(basque12$lr_pse, basque12$natid,mean)
means_lr_pse_natid09   <- tapply(basque09$lr_pse, basque09$natid,mean)

# put years together
means_lr_pse_natid <- cbind(round(means_lr_pse_natid20, digits=3), 
                            round(means_lr_pse_natid16, digits=3), 
                            round(means_lr_pse_natid12, digits=3),
                            round(means_lr_pse_natid09, digits=3))
means_lr_pse_natid

#---------
# ideo5
#---------
means_lr_pse_ideo520   <- tapply(basque20$lr_pse, basque20$ideo5,mean)
means_lr_pse_ideo516   <- tapply(basque16$lr_pse, basque16$ideo5,mean)
means_lr_pse_ideo512   <- tapply(basque12$lr_pse, basque12$ideo5,mean)
means_lr_pse_ideo509   <- tapply(basque09$lr_pse, basque09$ideo5,mean)

# put years together
means_lr_pse_ideo5 <- cbind(round(means_lr_pse_ideo520, digits=3),  
                            round(means_lr_pse_ideo516, digits=3),  
                            round(means_lr_pse_ideo512, digits=3),  
                            round(means_lr_pse_ideo509, digits=3))
means_lr_pse_ideo5

#---------
# gender
#---------
means_lr_pse_gender20   <- tapply(basque20$lr_pse, basque20$gender,mean)
means_lr_pse_gender16   <- tapply(basque16$lr_pse, basque16$gender,mean)
means_lr_pse_gender12   <- tapply(basque12$lr_pse, basque12$gender,mean)
means_lr_pse_gender09   <- tapply(basque09$lr_pse, basque09$gender,mean)


# put years together
means_lr_pse_gender <- cbind(round(means_lr_pse_gender20, digits=3),  
                             round(means_lr_pse_gender16, digits=3),  
                             round(means_lr_pse_gender12, digits=3),  
                             round(means_lr_pse_gender09, digits=3))
means_lr_pse_gender


#---------
# cath (not contained in 2009)
#---------
means_lr_pse_cath20   <- tapply(basque20$lr_pse, basque20$cath,mean)
means_lr_pse_cath16   <- tapply(basque16$lr_pse, basque16$cath,mean)
means_lr_pse_cath12   <- tapply(basque12$lr_pse, basque12$cath,mean)


# put years together
means_lr_pse_cath <- cbind(round(means_lr_pse_cath20, digits=3),  
                           round(means_lr_pse_cath16, digits=3),  
                           round(means_lr_pse_cath12, digits=3), NA)
means_lr_pse_cath

#---------
# status
#---------
means_lr_pse_status20   <- tapply(basque20$lr_pse, basque20$status,mean)
means_lr_pse_status16   <- tapply(basque16$lr_pse, basque16$status,mean)
means_lr_pse_status12   <- tapply(basque12$lr_pse, basque12$status,mean)
means_lr_pse_status09   <- tapply(basque09$lr_pse, basque09$status,mean)


# put years together
means_lr_pse_status <- cbind(round(means_lr_pse_status20, digits=3),  
                             round(means_lr_pse_status16, digits=3),  
                             round(means_lr_pse_status12, digits=3),  
                             round(means_lr_pse_status09, digits=3))
means_lr_pse_status

#---------------------
# put variables together
#---------------------
tableA3 <- rbind(means_lr_pse_natid,means_lr_pse_ideo5, 
                       means_lr_pse_gender, means_lr_pse_cath, means_lr_pse_status)
colnames(tableA3) <- c(2020, 2016, 2012, 2009)
tableA3



#------------------------------------------------------------------------------------
# Table A4: Mean Left-Right Party Placement Perceptions by Groups: PP 
#------------------------------------------------------------------------------------

#----------
# natid
#----------
means_lr_pp_natid20   <- tapply(basque20$lr_pp, basque20$natid,mean)
means_lr_pp_natid16   <- tapply(basque16$lr_pp, basque16$natid,mean)
means_lr_pp_natid12   <- tapply(basque12$lr_pp, basque12$natid,mean)
means_lr_pp_natid09   <- tapply(basque09$lr_pp, basque09$natid,mean)

# put years together
means_lr_pp_natid <- cbind(round(means_lr_pp_natid20, digits=3), 
                           round(means_lr_pp_natid16, digits=3), 
                           round(means_lr_pp_natid12, digits=3),
                           round(means_lr_pp_natid09, digits=3))
means_lr_pp_natid

#---------
# ideo5
#---------
means_lr_pp_ideo520   <- tapply(basque20$lr_pp, basque20$ideo5,mean)
means_lr_pp_ideo516   <- tapply(basque16$lr_pp, basque16$ideo5,mean)
means_lr_pp_ideo512   <- tapply(basque12$lr_pp, basque12$ideo5,mean)
means_lr_pp_ideo509   <- tapply(basque09$lr_pp, basque09$ideo5,mean)

# put years together
means_lr_pp_ideo5 <- cbind(round(means_lr_pp_ideo520, digits=3),  
                           round(means_lr_pp_ideo516, digits=3),  
                           round(means_lr_pp_ideo512, digits=3),  
                           round(means_lr_pp_ideo509, digits=3))
means_lr_pp_ideo5

#---------
# gender
#---------
means_lr_pp_gender20   <- tapply(basque20$lr_pp, basque20$gender,mean)
means_lr_pp_gender16   <- tapply(basque16$lr_pp, basque16$gender,mean)
means_lr_pp_gender12   <- tapply(basque12$lr_pp, basque12$gender,mean)
means_lr_pp_gender09   <- tapply(basque09$lr_pp, basque09$gender,mean)

# put years together
means_lr_pp_gender <- cbind(round(means_lr_pp_gender20, digits=3),  
                            round(means_lr_pp_gender16, digits=3),  
                            round(means_lr_pp_gender12, digits=3),  
                            round(means_lr_pp_gender09, digits=3))
means_lr_pp_gender

#---------
# cath (not contained in 2009)
#---------
means_lr_pp_cath20   <- tapply(basque20$lr_pp, basque20$cath,mean)
means_lr_pp_cath16   <- tapply(basque16$lr_pp, basque16$cath,mean)
means_lr_pp_cath12   <- tapply(basque12$lr_pp, basque12$cath,mean)

# put years together
means_lr_pp_cath <- cbind(round(means_lr_pp_cath20, digits=3),  
                          round(means_lr_pp_cath16, digits=3),  
                          round(means_lr_pp_cath12, digits=3), NA)
means_lr_pp_cath

#---------
# status
#---------
means_lr_pp_status20   <- tapply(basque20$lr_pp, basque20$status,mean)
means_lr_pp_status16   <- tapply(basque16$lr_pp, basque16$status,mean)
means_lr_pp_status12   <- tapply(basque12$lr_pp, basque12$status,mean)
means_lr_pp_status09   <- tapply(basque09$lr_pp, basque09$status,mean)

# put years together
means_lr_pp_status <- cbind(round(means_lr_pp_status20, digits=3),  
                            round(means_lr_pp_status16, digits=3),  
                            round(means_lr_pp_status12, digits=3),  
                            round(means_lr_pp_status09, digits=3))
means_lr_pp_status


#---------------------
# put variables together
#---------------------
tableA4 <- rbind(means_lr_pp_natid,means_lr_pp_ideo5, 
                      means_lr_pp_gender, means_lr_pp_cath, means_lr_pp_status)
colnames(tableA4) <- c(2020, 2016, 2012, 2009)
tableA4 


#------------------------------------------------------------------------------------
# Table A5: Mean Left-Right Party Placement Perceptions by Groups: EH-Bildu
#------------------------------------------------------------------------------------

#----------
# natid
#----------
means_lr_eh_natid20   <- tapply(basque20$lr_eh, basque20$natid,mean)
means_lr_eh_natid16   <- tapply(basque16$lr_eh, basque16$natid,mean)
means_lr_eh_natid12   <- tapply(basque12$lr_eh, basque12$natid,mean)
means_lr_eh_natid09   <- tapply(basque09$lr_aralar, basque09$natid,mean)

# put years together
means_lr_eh_natid <- cbind(round(means_lr_eh_natid20, digits=3), 
                           round(means_lr_eh_natid16, digits=3), 
                           round(means_lr_eh_natid12, digits=3),
                           round(means_lr_eh_natid09, digits=3))
means_lr_eh_natid

#---------
# ideo5
#---------
means_lr_eh_ideo520   <- tapply(basque20$lr_eh, basque20$ideo5,mean)
means_lr_eh_ideo516   <- tapply(basque16$lr_eh, basque16$ideo5,mean)
means_lr_eh_ideo512   <- tapply(basque12$lr_eh, basque12$ideo5,mean)
means_lr_eh_ideo509   <- tapply(basque09$lr_aralar, basque09$ideo5,mean)

# put years together
means_lr_eh_ideo5 <- cbind(round(means_lr_eh_ideo520, digits=3),  
                           round(means_lr_eh_ideo516, digits=3),  
                           round(means_lr_eh_ideo512, digits=3),  
                           round(means_lr_eh_ideo509, digits=3))
means_lr_eh_ideo5

#---------
# gender
#---------
means_lr_eh_gender20   <- tapply(basque20$lr_eh, basque20$gender,mean)
means_lr_eh_gender16   <- tapply(basque16$lr_eh, basque16$gender,mean)
means_lr_eh_gender12   <- tapply(basque12$lr_eh, basque12$gender,mean)
means_lr_eh_gender09   <- tapply(basque09$lr_aralar, basque09$gender,mean)

# put years together
means_lr_eh_gender <- cbind(round(means_lr_eh_gender20, digits=3),  
                            round(means_lr_eh_gender16, digits=3),  
                            round(means_lr_eh_gender12, digits=3),  
                            round(means_lr_eh_gender09, digits=3))
means_lr_eh_gender

#---------
# cath (not contained in 2009)
#---------
means_lr_eh_cath20   <- tapply(basque20$lr_eh, basque20$cath,mean)
means_lr_eh_cath16   <- tapply(basque16$lr_eh, basque16$cath,mean)
means_lr_eh_cath12   <- tapply(basque12$lr_eh, basque12$cath,mean)

# put years together
means_lr_eh_cath <- cbind(round(means_lr_eh_cath20, digits=3),  
                          round(means_lr_eh_cath16, digits=3),  
                          round(means_lr_eh_cath12, digits=3), NA)
means_lr_eh_cath


#---------
# status
#---------
means_lr_eh_status20   <- tapply(basque20$lr_eh, basque20$status,mean)
means_lr_eh_status16   <- tapply(basque16$lr_eh, basque16$status,mean)
means_lr_eh_status12   <- tapply(basque12$lr_eh, basque12$status,mean)
means_lr_eh_status09   <- tapply(basque09$lr_aralar, basque09$status,mean)

# put years together
means_lr_eh_status <- cbind(round(means_lr_eh_status20, digits=3),  
                            round(means_lr_eh_status16, digits=3),  
                            round(means_lr_eh_status12, digits=3),  
                            round(means_lr_eh_status09, digits=3))
means_lr_eh_status


#---------------------
# put variables together
#---------------------
tableA5 <- rbind(means_lr_eh_natid,means_lr_eh_ideo5, 
                      means_lr_eh_gender, means_lr_eh_cath, means_lr_eh_status)
colnames(tableA5) <- c(2020, 2016, 2012, 2009)
tableA5


#------------------------------------------------------------------------------------
# Table A6: Mean Left-Right Party Placement Perceptions by Groups: PNV 
#------------------------------------------------------------------------------------

#----------
# natid
#----------
means_lr_pnv_natid20   <- tapply(basque20$lr_pnv, basque20$natid,mean)
means_lr_pnv_natid16   <- tapply(basque16$lr_pnv, basque16$natid,mean)
means_lr_pnv_natid12   <- tapply(basque12$lr_pnv, basque12$natid,mean)
means_lr_pnv_natid09   <- tapply(basque09$lr_pnv, basque09$natid,mean)

# put years together
means_lr_pnv_natid <- cbind(round(means_lr_pnv_natid20, digits=3), 
                            round(means_lr_pnv_natid16, digits=3), 
                            round(means_lr_pnv_natid12, digits=3),
                            round(means_lr_pnv_natid09, digits=3))
means_lr_pnv_natid

#---------
# ideo5
#---------
means_lr_pnv_ideo520   <- tapply(basque20$lr_pnv, basque20$ideo5,mean)
means_lr_pnv_ideo516   <- tapply(basque16$lr_pnv, basque16$ideo5,mean)
means_lr_pnv_ideo512   <- tapply(basque12$lr_pnv, basque12$ideo5,mean)
means_lr_pnv_ideo509   <- tapply(basque09$lr_pnv, basque09$ideo5,mean)

# put years together
means_lr_pnv_ideo5 <- cbind(round(means_lr_pnv_ideo520, digits=3),  
                            round(means_lr_pnv_ideo516, digits=3),  
                            round(means_lr_pnv_ideo512, digits=3),  
                            round(means_lr_pnv_ideo509, digits=3))
means_lr_pnv_ideo5

#---------
# gender
#---------
means_lr_pnv_gender20   <- tapply(basque20$lr_pnv, basque20$gender,mean)
means_lr_pnv_gender16   <- tapply(basque16$lr_pnv, basque16$gender,mean)
means_lr_pnv_gender12   <- tapply(basque12$lr_pnv, basque12$gender,mean)
means_lr_pnv_gender09   <- tapply(basque09$lr_pnv, basque09$gender,mean)

# put years together
means_lr_pnv_gender <- cbind(round(means_lr_pnv_gender20, digits=3),  
                             round(means_lr_pnv_gender16, digits=3),  
                             round(means_lr_pnv_gender12, digits=3),  
                             round(means_lr_pnv_gender09, digits=3))
means_lr_pnv_gender

#---------
# cath (not contained in 2009)
#---------
means_lr_pnv_cath20   <- tapply(basque20$lr_pnv, basque20$cath,mean)
means_lr_pnv_cath16   <- tapply(basque16$lr_pnv, basque16$cath,mean)
means_lr_pnv_cath12   <- tapply(basque12$lr_pnv, basque12$cath,mean)

# put years together
means_lr_pnv_cath <- cbind(round(means_lr_pnv_cath20, digits=3),  
                           round(means_lr_pnv_cath16, digits=3),  
                           round(means_lr_pnv_cath12, digits=3), NA)
means_lr_pnv_cath

#---------
# status
#---------
means_lr_pnv_status20   <- tapply(basque20$lr_pnv, basque20$status,mean)
means_lr_pnv_status16   <- tapply(basque16$lr_pnv, basque16$status,mean)
means_lr_pnv_status12   <- tapply(basque12$lr_pnv, basque12$status,mean)
means_lr_pnv_status09   <- tapply(basque09$lr_pnv, basque09$status,mean)

# put years together
means_lr_pnv_status <- cbind(round(means_lr_pnv_status20, digits=3),  
                             round(means_lr_pnv_status16, digits=3),  
                             round(means_lr_pnv_status12, digits=3),  
                             round(means_lr_pnv_status09, digits=3))
means_lr_pnv_status


#---------------------
# put variables together
#---------------------
tableA6 <- rbind(means_lr_pnv_natid,means_lr_pnv_ideo5, 
                       means_lr_pnv_gender, means_lr_pnv_cath, means_lr_pnv_status)
colnames(tableA6) <- c(2020, 2016, 2012, 2009)
tableA6



#===========================================
# D Estimation Tables and Coefficient Plots
#===========================================


#------------------------------------------------------------------------------------
# Table A7: Multivariate Linear Model Estimates, 2020 Basque Regional Election
#------------------------------------------------------------------------------------

#---------------------------
# save coefs. in a matrix
#---------------------------
coef.m20<-coef(m20)
colnames(coef.m20)<-c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
rownames(coef.m20)<-c("Intercept", "National Sentiment: Spanish", "National Sentiment: Basque", "Basque Language Usage",
                      "Left-Right Self-placement", 
                      "Ideological Self-categorization: Conservative", "Ideological Self-categorization: Socialist",
                      "Ideological Self-categorization: Feminist.Ecologist", "Ideological Self-categorization: Liberal",
                      "Economic Situation","Female","Age", "Social Class: high", "Social Class: low", "Catholic", "Urbanization")
coef.m20


#----------------------------------
# save standard errors in a matrix
#----------------------------------
set_m20<- rbind(summary(m20)$`Response lr_podemos`$coefficients[,2],
                summary(m20)$`Response lr_pse`$coefficients[,2],
                summary(m20)$`Response lr_pp`$coefficients[,2], 
                summary(m20)$`Response lr_eh`$coefficients[,2],
                summary(m20)$`Response lr_pnv`$coefficients[,2])

se_m20 <- t(set_m20)
colnames(se_m20)<-c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
rownames(se_m20)<-c("Intercept", "National Sentiment: Spanish", "National Sentiment: Basque", "Basque Language Usage",
                    "Left-Right Self-placement", 
                    "Ideological Self-categorization: Conservative", "Ideological Self-categorization: Socialist",
                    "Ideological Self-categorization: Feminist.Ecologist", "Ideological Self-categorization: Liberal",
                    "Economic Situation","Female","Age", "Social Class: high", "Social Class: low", "Catholic", "Urbanization")

#------------------------------------
# save R squared values in a matrix
#------------------------------------
R_m20<-matrix(NA,1,5)
colnames(R_m20)<-c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
rownames(R_m20)<-c("r.squared")

R_m20["r.squared", "Podemos"]  <-summary(m20)$`Response lr_podemos`$r.squared
R_m20["r.squared", "PSE"]      <-summary(m20)$`Response lr_pse`$r.squared
R_m20["r.squared", "PP"]       <-summary(m20)$`Response lr_pp`$r.squared
R_m20["r.squared", "EH-Bildu"] <-summary(m20)$`Response lr_eh`$r.squared
R_m20["r.squared", "PNV"]      <-summary(m20)$`Response lr_pnv`$r.squared

#---------------------------
# put the matrices together
#---------------------------
coef.se.p_m20<-matrix(NA,5,32)
rownames(coef.se.p_m20)<-c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
colnames(coef.se.p_m20)<-c("National Sentiment: Spanish", "s.e.",
                               "National Sentiment: Basque","s.e.",  
                               "Basque Language Usage","s.e.", 
                               "Left-Right Self-placement","s.e.",  
                               "Ideological Self-categorization: Conservative", "s.e.",
                               "Ideological Self-categorization: Socialist","s.e.",
                               "Ideological Self-categorization: Feminist.Ecologist", "s.e.",
                               "Ideological Self-categorization: Liberal","s.e.",
                               "Economic Situation","s.e.", 
                               "Female","s.e.", 
                               "Age", "s.e.",
                               "Social Class: high", "s.e.", 
                               "Social Class: low", "s.e.", 
                               "Catholic", "s.e.", 
                               "Urbanization", "s.e.",
                               "Intercept","s.e.")


# natidSpanish:
coef.se.p_m20[,1]  <-t(coef.m20[2,1:5])
coef.se.p_m20[,2]  <-t(se_m20[2,1:5])

# natidBasque:
coef.se.p_m20[,3]  <-t(coef.m20[3,1:5])
coef.se.p_m20[,4]  <-t(se_m20[3,1:5])

# eusk:
coef.se.p_m20[,5]  <-t(coef.m20[4,1:5])
coef.se.p_m20[,6]  <-t(se_m20[4,1:5])

# lr_ego:
coef.se.p_m20[,7]  <-t(coef.m20[5,1:5])
coef.se.p_m20[,8]  <-t(se_m20[5,1:5])

# ideo5Conservative:
coef.se.p_m20[,9]  <-t(coef.m20[6,1:5])
coef.se.p_m20[,10] <-t(se_m20[6,1:5])

# ideo5Socialist:
coef.se.p_m20[,11] <-t(coef.m20[7,1:5])
coef.se.p_m20[,12] <-t(se_m20[7,1:5])

# ideo5Feminist.Ecologist:
coef.se.p_m20[,13] <-t(coef.m20[8,1:5])
coef.se.p_m20[,14] <-t(se_m20[8,1:5])

# ideo5Liberal:
coef.se.p_m20[,15] <-t(coef.m20[9,1:5])
coef.se.p_m20[,16] <-t(se_m20[9,1:5])

# eco:
coef.se.p_m20[,17] <-t(coef.m20[10,1:5])
coef.se.p_m20[,18] <-t(se_m20[10,1:5])

# gender:
coef.se.p_m20[,19] <-t(coef.m20[11,1:5])
coef.se.p_m20[,20] <-t(se_m20[11,1:5])

# age:
coef.se.p_m20[,21] <-t(coef.m20[12,1:5])
coef.se.p_m20[,22] <-t(se_m20[12,1:5])

# statuslow:
coef.se.p_m20[,23] <-t(coef.m20[13,1:5])
coef.se.p_m20[,24] <-t(se_m20[13,1:5])

# statushigh:
coef.se.p_m20[,25] <-t(coef.m20[14,1:5])
coef.se.p_m20[,26] <-t(se_m20[14,1:5])

# catholic
coef.se.p_m20[,27] <-t(coef.m20[15,1:5])
coef.se.p_m20[,28] <-t(se_m20[15,1:5])

# pop
coef.se.p_m20[,29] <-t(coef.m20[16,1:5])
coef.se.p_m20[,30] <-t(se_m20[16,1:5])

# intercept: 
coef.se.p_m20[,31] <-t(coef.m20[1,1:5])
coef.se.p_m20[,32] <-t(se_m20[1,1:5])

coef.se.p_m20

# transpose matrix
Tcoef.se.p_m20<-t(coef.se.p_m20)

fitTm20<-matrix(NA,1,5)
colnames(fitTm20)<-c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
rownames(fitTm20)<-c("Model Fit")
fitTm20

tableA7 <-rbind(Tcoef.se.p_m20, fitTm20, R_m20)
tableA7


#------------------------------------------------------------------------------------
# Table A8: Multivariate Linear Model Estimates, 2016 Basque Regional Election
#------------------------------------------------------------------------------------

#---------------------------
# save coefs. in a matrix
#---------------------------
coef.m16<-coef(m16)
colnames(coef.m16) <-c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
rownames(coef.m16) <-c("Intercept", "National Sentiment: Spanish", "National Sentiment: Basque", "Basque Language Usage",
                        "Left-Right Self-placement", 
                        "Ideological Self-categorization: Conservative", "Ideological Self-categorization: Socialist",
                        "Ideological Self-categorization: Feminist.Ecologist", "Ideological Self-categorization: Liberal",
                        "Economic Situation","Female","Age", "Social Class: high", "Social Class: low", "Catholic", "Urbanization")

#----------------------------------
# save standard errors in a matrix
#----------------------------------
set_m16<- rbind(summary(m16)$`Response lr_podemos`$coefficients[,2],
                summary(m16)$`Response lr_pse`$coefficients[,2],
                summary(m16)$`Response lr_pp`$coefficients[,2], 
                summary(m16)$`Response lr_eh`$coefficients[,2],
                summary(m16)$`Response lr_pnv`$coefficients[,2])

se_m16 <- t(set_m16)
colnames(se_m16) <-c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
rownames(se_m16) <-c("Intercept", "National Sentiment: Spanish", "National Sentiment: Basque", "Basque Language Usage",
                      "Left-Right Self-placement", 
                      "Ideological Self-categorization: Conservative", "Ideological Self-categorization: Socialist",
                      "Ideological Self-categorization: Feminist.Ecologist", "Ideological Self-categorization: Liberal",
                      "Economic Situation","Female","Age", "Social Class: high", "Social Class: low", "Catholic", "Urbanization")

#------------------------------------
# save R squared values in a matrix
#------------------------------------
R_m16<-matrix(NA,1,5)
colnames(R_m16) <-c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
rownames(R_m16) <-c("r.squared")

R_m16["r.squared", "Podemos"]   <-summary(m16)$`Response lr_podemos`$r.squared
R_m16["r.squared", "PSE"]       <-summary(m16)$`Response lr_pse`$r.squared
R_m16["r.squared", "PP"]        <-summary(m16)$`Response lr_pp`$r.squared
R_m16["r.squared", "EH-Bildu"]  <-summary(m16)$`Response lr_eh`$r.squared
R_m16["r.squared", "PNV"]       <-summary(m16)$`Response lr_pnv`$r.squared

#---------------------------
# put the matrices together
#---------------------------
coef.se.p_m16<-matrix(NA,5,32)
rownames(coef.se.p_m16)<-c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
colnames(coef.se.p_m16)<-c("National Sentiment: Spanish", "s.e.",
                            "National Sentiment: Basque","s.e.",  
                            "Basque Language Usage","s.e.", 
                            "Left-Right Self-placement","s.e.",  
                            "Ideological Self-categorization: Conservative", "s.e.",
                            "Ideological Self-categorization: Socialist","s.e.",
                            "Ideological Self-categorization: Feminist.Ecologist", "s.e.",
                            "Ideological Self-categorization: Liberal","s.e.",
                            "Economic Situation","s.e.", 
                            "Female","s.e.", 
                            "Age", "s.e.",
                            "Social Class: high", "s.e.", 
                            "Social Class: low", "s.e.", 
                            "Catholic", "s.e.", 
                            "Urbanization", "s.e.",
                            "Intercept","s.e.")


# natidSpanish:
coef.se.p_m16[,1] <-t(coef.m16[2,1:5])
coef.se.p_m16[,2] <-t(se_m16[2,1:5])

# natidBasque:
coef.se.p_m16[,3] <-t(coef.m16[3,1:5])
coef.se.p_m16[,4] <-t(se_m16[3,1:5])

# eusk:
coef.se.p_m16[,5] <-t(coef.m16[4,1:5])
coef.se.p_m16[,6] <-t(se_m16[4,1:5])

# lr_ego:
coef.se.p_m16[,7] <-t(coef.m16[5,1:5])
coef.se.p_m16[,8] <-t(se_m16[5,1:5])

# ideo5Conservative:
coef.se.p_m16[,9] <-t(coef.m16[6,1:5])
coef.se.p_m16[,10] <-t(se_m16[6,1:5])

# ideo5Socialist:
coef.se.p_m16[,11] <-t(coef.m16[7,1:5])
coef.se.p_m16[,12] <-t(se_m16[7,1:5])

# ideo5Feminist.Ecologist:
coef.se.p_m16[,13] <-t(coef.m16[8,1:5])
coef.se.p_m16[,14] <-t(se_m16[8,1:5])

# ideo5Liberal:
coef.se.p_m16[,15] <-t(coef.m16[9,1:5])
coef.se.p_m16[,16] <-t(se_m16[9,1:5])

# eco:
coef.se.p_m16[,17] <-t(coef.m16[10,1:5])
coef.se.p_m16[,18] <-t(se_m16[10,1:5])

# gender:
coef.se.p_m16[,19] <-t(coef.m16[11,1:5])
coef.se.p_m16[,20] <-t(se_m16[11,1:5])

# age:
coef.se.p_m16[,21] <-t(coef.m16[12,1:5])
coef.se.p_m16[,22] <-t(se_m16[12,1:5])

# statuslow:
coef.se.p_m16[,23] <-t(coef.m16[13,1:5])
coef.se.p_m16[,24] <-t(se_m16[13,1:5])

# statushigh:
coef.se.p_m16[,25] <-t(coef.m16[14,1:5])
coef.se.p_m16[,26] <-t(se_m16[14,1:5])

# catholic
coef.se.p_m16[,27] <-t(coef.m16[15,1:5])
coef.se.p_m16[,28] <-t(se_m16[15,1:5])

# pop
coef.se.p_m16[,29] <-t(coef.m16[16,1:5])
coef.se.p_m16[,30] <-t(se_m16[16,1:5])

# intercept: 
coef.se.p_m16[,31] <-t(coef.m16[1,1:5])
coef.se.p_m16[,32] <-t(se_m16[1,1:5])


# transpose matrix
Tcoef.se.p_m16<-t(coef.se.p_m16)

fitTm16<-matrix(NA,1,5)
colnames(fitTm16)<-c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
rownames(fitTm16)<-c("Model Fit")

tableA8 <-rbind(Tcoef.se.p_m16, fitTm16, R_m16)
tableA8



#------------------------------------------------------------------------------------
# Table A9: Multivariate Linear Model Estimates, 2012 Basque Regional Election 
#------------------------------------------------------------------------------------

#---------------------------
# save coefs. in a matrix
#---------------------------
coef.m12<-coef(m12)
colnames(coef.m12)<-c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
rownames(coef.m12)<-c("Intercept", "National Sentiment: Spanish", "National Sentiment: Basque", "Basque Language Usage",
                      "Left-Right Self-placement", 
                      "Ideological Self-categorization: Conservative", "Ideological Self-categorization: Socialist",
                      "Ideological Self-categorization: Feminist.Ecologist", "Ideological Self-categorization: Liberal",
                      "Economic Situation","Female","Age", "Social Class: high", "Social Class: low", "Catholic", "Urbanization")

#----------------------------------
# save standard errors in a matrix
#----------------------------------
set_m12<- rbind(summary(m12)$`Response lr_iu`$coefficients[,2],
                summary(m12)$`Response lr_pse`$coefficients[,2],
                summary(m12)$`Response lr_pp`$coefficients[,2], 
                summary(m12)$`Response lr_eh`$coefficients[,2],
                summary(m12)$`Response lr_pnv`$coefficients[,2])

se_m12 <- t(set_m12)
colnames(se_m12)<-c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
rownames(se_m12)<-c("Intercept", "National Sentiment: Spanish", "National Sentiment: Basque", "Basque Language Usage",
                    "Left-Right Self-placement", 
                    "Ideological Self-categorization: Conservative", "Ideological Self-categorization: Socialist",
                    "Ideological Self-categorization: Feminist.Ecologist", "Ideological Self-categorization: Liberal",
                    "Economic Situation","Female","Age", "Social Class: high", "Social Class: low", "Catholic", "Urbanization")

#------------------------------------
# save R squared values in a matrix
#------------------------------------
R_m12<-matrix(NA,1,5)
colnames(R_m12)<-c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
rownames(R_m12)<-c("r.squared")

R_m12["r.squared", "Podemos"]  <-summary(m12)$`Response lr_iu`$r.squared
R_m12["r.squared", "PSE"]      <-summary(m12)$`Response lr_pse`$r.squared
R_m12["r.squared", "PP"]       <-summary(m12)$`Response lr_pp`$r.squared
R_m12["r.squared", "EH-Bildu"] <-summary(m12)$`Response lr_eh`$r.squared
R_m12["r.squared", "PNV"]      <-summary(m12)$`Response lr_pnv`$r.squared

#---------------------------
# put the matrices together
#---------------------------
coef.se.p_m12<-matrix(NA,5,32)
rownames(coef.se.p_m12)<-c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
colnames(coef.se.p_m12)<-c("National Sentiment: Spanish", "s.e.",
                            "National Sentiment: Basque","s.e.",  
                            "Basque Language Usage","s.e.", 
                            "Left-Right Self-placement","s.e.",  
                            "Ideological Self-categorization: Conservative", "s.e.",
                            "Ideological Self-categorization: Socialist","s.e.",
                            "Ideological Self-categorization: Feminist.Ecologist", "s.e.",
                            "Ideological Self-categorization: Liberal","s.e.",
                            "Economic Situation","s.e.", 
                            "Female","s.e.", 
                            "Age", "s.e.",
                            "Social Class: high", "s.e.", 
                            "Social Class: low", "s.e.", 
                            "Catholic", "s.e.", 
                            "Urbanization", "s.e.",
                            "Intercept","s.e.")


# natidSpanish:
coef.se.p_m12[,1] <-t(coef.m12[2,1:5])
coef.se.p_m12[,2] <-t(se_m12[2,1:5])

# natidBasque:
coef.se.p_m12[,3] <-t(coef.m12[3,1:5])
coef.se.p_m12[,4] <-t(se_m12[3,1:5])

# eusk:
coef.se.p_m12[,5] <-t(coef.m12[4,1:5])
coef.se.p_m12[,6] <-t(se_m12[4,1:5])

# lr_ego:
coef.se.p_m12[,7] <-t(coef.m12[5,1:5])
coef.se.p_m12[,8] <-t(se_m12[5,1:5])

# ideo5Conservative:
coef.se.p_m12[,9] <-t(coef.m12[6,1:5])
coef.se.p_m12[,10] <-t(se_m12[6,1:5])

# ideo5Socialist:
coef.se.p_m12[,11] <-t(coef.m12[7,1:5])
coef.se.p_m12[,12] <-t(se_m12[7,1:5])

# ideo5Feminist.Ecologist:
coef.se.p_m12[,13] <-t(coef.m12[8,1:5])
coef.se.p_m12[,14] <-t(se_m12[8,1:5])

# ideo5Liberal:
coef.se.p_m12[,15] <-t(coef.m12[9,1:5])
coef.se.p_m12[,16] <-t(se_m12[9,1:5])

# eco:
coef.se.p_m12[,17] <-t(coef.m12[10,1:5])
coef.se.p_m12[,18] <-t(se_m12[10,1:5])

# gender:
coef.se.p_m12[,19] <-t(coef.m12[11,1:5])
coef.se.p_m12[,20] <-t(se_m12[11,1:5])

# age:
coef.se.p_m12[,21] <-t(coef.m12[12,1:5])
coef.se.p_m12[,22] <-t(se_m12[12,1:5])

# statuslow:
coef.se.p_m12[,23] <-t(coef.m12[13,1:5])
coef.se.p_m12[,24] <-t(se_m12[13,1:5])

# statushigh:
coef.se.p_m12[,25] <-t(coef.m12[14,1:5])
coef.se.p_m12[,26] <-t(se_m12[14,1:5])

# catholic
coef.se.p_m12[,27] <-t(coef.m12[15,1:5])
coef.se.p_m12[,28] <-t(se_m12[15,1:5])

# pop
coef.se.p_m12[,29] <-t(coef.m12[16,1:5])
coef.se.p_m12[,30] <-t(se_m12[16,1:5])

# intercept: 
coef.se.p_m12[,31] <-t(coef.m12[1,1:5])
coef.se.p_m12[,32] <-t(se_m12[1,1:5])

coef.se.p_m12

# transpose matrix
Tcoef.se.p_m12<-t(coef.se.p_m12)

fitTm12<-matrix(NA,1,5)
colnames(fitTm12)<-c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
rownames(fitTm12)<-c("Model Fit")

tableA9 <-rbind(Tcoef.se.p_m12, fitTm12, R_m12)
tableA9


#------------------------------------------------------------------------------------
# Table A10: Multivariate Linear Model Estimates, 2009 Basque Regional Election 
#------------------------------------------------------------------------------------

#---------------------------
# save coefs. in a matrix
#---------------------------
coef.m09<-coef(m09)
colnames(coef.m09) <-c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
rownames(coef.m09) <-c("Intercept", "National Sentiment: Spanish", "National Sentiment: Basque", "Basque Language Usage",
                        "Left-Right Self-placement", 
                        "Ideological Self-categorization: Conservative", "Ideological Self-categorization: Socialist",
                        "Ideological Self-categorization: Feminist.Ecologist", "Ideological Self-categorization: Liberal",
                        "Economic Situation","Female","Age", "Social Class: high", "Social Class: low", "Urbanization")

#----------------------------------
# save standard errors in a matrix
#----------------------------------
set_m09<- rbind(summary(m09)$`Response lr_eb`$coefficients[,2],
                summary(m09)$`Response lr_pse`$coefficients[,2],
                summary(m09)$`Response lr_pp`$coefficients[,2], 
                summary(m09)$`Response lr_aralar`$coefficients[,2],
                summary(m09)$`Response lr_pnv`$coefficients[,2])

se_m09 <- t(set_m09)
colnames(se_m09)<-c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
rownames(se_m09)<-c("Intercept", "National Sentiment: Spanish", "National Sentiment: Basque", "Basque Language Usage",
                    "Left-Right Self-placement", 
                    "Ideological Self-categorization: Conservative", "Ideological Self-categorization: Socialist",
                    "Ideological Self-categorization: Feminist.Ecologist", "Ideological Self-categorization: Liberal",
                    "Economic Situation","Female","Age", "Social Class: high", "Social Class: low", "Urbanization")


#------------------------------------
# save R squared values in a matrix
#------------------------------------
R_m09<-matrix(NA,1,5)
colnames(R_m09)<-c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
rownames(R_m09)<-c("r.squared")

R_m09["r.squared","Podemos"]   <-summary(m09)$`Response lr_eb`$r.squared
R_m09["r.squared","PSE"]       <-summary(m09)$`Response lr_pse`$r.squared
R_m09["r.squared","PP"]        <-summary(m09)$`Response lr_pp`$r.squared
R_m09["r.squared","EH-Bildu"]  <-summary(m09)$`Response lr_aralar`$r.squared
R_m09["r.squared","PNV"]       <-summary(m09)$`Response lr_pnv`$r.squared

#---------------------------
# put the matrices together
#---------------------------

coef.se.p_m09<-matrix(NA,5,30)
rownames(coef.se.p_m09)<-c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
colnames(coef.se.p_m09)<-c("National Sentiment: Spanish", "s.e.",
                            "National Sentiment: Basque","s.e.",  
                            "Basque Language Usage","s.e.", 
                            "Left-Right Self-placement","s.e.",  
                            "Ideological Self-categorization: Conservative", "s.e.",
                            "Ideological Self-categorization: Socialist","s.e.",
                            "Ideological Self-categorization: Feminist.Ecologist", "s.e.",
                            "Ideological Self-categorization: Liberal","s.e.",
                            "Economic Situation","s.e.", 
                            "Female","s.e.", 
                            "Age", "s.e.",
                            "Social Class: high", "s.e.", 
                            "Social Class: low", "s.e.", 
                            "Urbanization", "s.e.",
                            "Intercept","s.e.")


# natidSpanish:
coef.se.p_m09[,1] <-t(coef.m09[2,1:5])
coef.se.p_m09[,2] <-t(se_m09[2,1:5])

# natidBasque:
coef.se.p_m09[,3] <-t(coef.m09[3,1:5])
coef.se.p_m09[,4] <-t(se_m09[3,1:5])

# eusk:
coef.se.p_m09[,5] <-t(coef.m09[4,1:5])
coef.se.p_m09[,6] <-t(se_m09[4,1:5])

# lr_ego:
coef.se.p_m09[,7] <-t(coef.m09[5,1:5])
coef.se.p_m09[,8] <-t(se_m09[5,1:5])

# ideo5Conservative:
coef.se.p_m09[,9] <-t(coef.m09[6,1:5])
coef.se.p_m09[,10] <-t(se_m09[6,1:5])

# ideo5Socialist:
coef.se.p_m09[,11] <-t(coef.m09[7,1:5])
coef.se.p_m09[,12] <-t(se_m09[7,1:5])

# ideo5Feminist.Ecologist:
coef.se.p_m09[,13] <-t(coef.m09[8,1:5])
coef.se.p_m09[,14] <-t(se_m09[8,1:5])

# ideo5Liberal:
coef.se.p_m09[,15] <-t(coef.m09[9,1:5])
coef.se.p_m09[,16] <-t(se_m09[9,1:5])

# eco:
coef.se.p_m09[,17] <-t(coef.m09[10,1:5])
coef.se.p_m09[,18] <-t(se_m09[10,1:5])

# gender:
coef.se.p_m09[,19] <-t(coef.m09[11,1:5])
coef.se.p_m09[,20] <-t(se_m09[11,1:5])

# age:
coef.se.p_m09[,21] <-t(coef.m09[12,1:5])
coef.se.p_m09[,22] <-t(se_m09[12,1:5])

# statuslow:
coef.se.p_m09[,23] <-t(coef.m09[13,1:5])
coef.se.p_m09[,24] <-t(se_m09[13,1:5])

# statushigh:
coef.se.p_m09[,25] <-t(coef.m09[14,1:5])
coef.se.p_m09[,26] <-t(se_m09[14,1:5])

# pop
coef.se.p_m09[,27] <-t(coef.m09[15,1:5])
coef.se.p_m09[,28] <-t(se_m09[15,1:5])

# intercept: 
coef.se.p_m09[,29] <-t(coef.m09[1,1:5])
coef.se.p_m09[,30] <-t(se_m09[1,1:5])

coef.se.p_m09

# transpose matrix
Tcoef.se.p_m09<-t(coef.se.p_m09)

fitTm09<-matrix(NA,1,5)
colnames(fitTm09)<-c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
rownames(fitTm09)<-c("Model Fit")

tableA10 <-rbind(Tcoef.se.p_m09, fitTm09, R_m09)
tableA10


#=================================================
# Figure A1: Model Estimates on National Groups
#=================================================

# 95% confidence intervals
interval <- -qnorm((1-0.95)/2)

require(gridExtra)

# 2020

mat_natid_m20 <- matrix(NA, 10, ncol = 4)
colnames(mat_natid_m20) <- c("Variable","Coefficient", "SE", "NatID")
mat_natid_m20[,1] <- c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
mat_natid_m20[,2] <- t(coef.m20[2:3,1:5])
mat_natid_m20[,3] <- t(se_m20[2:3,1:5]) 
mat_natid_m20[,4] <- rep(c("Spanish", "Basque"), each=5) 
mat_natid_m20

natid_m20<-as.data.frame(mat_natid_m20)
natid_m20$Coefficient <-as.numeric(natid_m20$Coefficient)
natid_m20$SE <-as.numeric(natid_m20$SE)
natid_m20
natid_m20$Variable <- factor(natid_m20$Variable, c("PNV", "EH-Bildu", "PP", "PSE",  "Podemos"))


# 2016

mat_natid_m16 <- matrix(NA, 10, ncol = 4)
colnames(mat_natid_m16) <- c("Variable","Coefficient", "SE", "NatID")
mat_natid_m16[,1] <- c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
mat_natid_m16[,2] <- t(coef.m16[2:3,1:5])
mat_natid_m16[,3] <- t(se_m16[2:3,1:5]) 
mat_natid_m16[,4] <- rep(c("Spanish", "Basque"), each=5) 
mat_natid_m16

natid_m16<-as.data.frame(mat_natid_m16)
natid_m16$Coefficient <-as.numeric(natid_m16$Coefficient)
natid_m16$SE <-as.numeric(natid_m16$SE)
natid_m16
natid_m16$Variable <- factor(natid_m16$Variable, c("PNV", "EH-Bildu", "PP", "PSE",  "Podemos"))


# 2012:

mat_natid_m12 <- matrix(NA, 10, ncol = 4)
colnames(mat_natid_m12) <- c("Variable","Coefficient", "SE", "NatID")
mat_natid_m12[,1] <- c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
mat_natid_m12[,2] <- t(coef.m12[2:3,1:5])
mat_natid_m12[,3] <- t(se_m12[2:3,1:5]) 
mat_natid_m12[,4] <- rep(c("Spanish", "Basque"), each=5) 
mat_natid_m12

natid_m12<-as.data.frame(mat_natid_m12)
natid_m12$Coefficient <-as.numeric(natid_m12$Coefficient)
natid_m12$SE <-as.numeric(natid_m12$SE)
natid_m12
natid_m12$Variable <- factor(natid_m12$Variable, c("PNV", "EH-Bildu", "PP", "PSE",  "Podemos"))


# 2009

mat_natid_m09 <- matrix(NA, 10, ncol = 4)
colnames(mat_natid_m09) <- c("Variable","Coefficient", "SE", "NatID")
mat_natid_m09[,1] <- c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
mat_natid_m09[,2] <- t(coef.m09[2:3,1:5])
mat_natid_m09[,3] <- t(se_m09[2:3,1:5]) 
mat_natid_m09[,4] <- rep(c("Spanish", "Basque"), each=5) 
mat_natid_m09

natid_m09<-as.data.frame(mat_natid_m09)
natid_m09$Coefficient <-as.numeric(natid_m09$Coefficient)
natid_m09$SE <-as.numeric(natid_m09$SE)
natid_m09
natid_m09$Variable <- factor(natid_m09$Variable, c("PNV", "EH-Bildu", "PP", "PSE",  "Podemos"))


#---------------------
# create plots
#---------------------
# 2020: 
p_natid_m20<- ggplot(natid_m20, aes(color=NatID, group=NatID)) +  ggtitle("2020") +
  geom_pointrange(aes(x = Variable, y = Coefficient,  group=NatID, 
                      ymin = Coefficient - SE*interval, ymax = Coefficient + SE*interval),
                  position = position_dodge(width = 1/1.6), shape=rep(c(15,17), each=5)) +
  coord_flip() + theme_classic() + labs(y="", x="") +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_vline(xintercept = 2.5, color = gray(1/2), lty = 2) +
  guides(color=guide_legend(title="(ref: neutral)")) +
  theme(plot.title = element_text(size = 10),
        legend.position = "none",
        legend.text=element_text(size=9),
        legend.title=element_text(size=9 )) +
  scale_color_manual(values=c("grey47", "black")) +
  scale_y_continuous(name="", limits=c(-1.2,1)) 



# 2016
p_natid_m16<- ggplot(natid_m16, aes(color=NatID, group=NatID)) +  ggtitle("2016") +
  geom_pointrange(aes(x = Variable, y = Coefficient,  group=NatID, 
                      ymin = Coefficient - SE*interval, ymax = Coefficient + SE*interval),
                  position = position_dodge(width = 1/1.6), shape=rep(c(15,17), each=5)) +
  coord_flip() + theme_classic() + labs(y="", x="") +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_vline(xintercept = 2.5, color = gray(1/2), lty = 2) +
  guides(color=guide_legend(title="(ref: neutral)")) +
  theme(plot.title = element_text(size = 10),
        legend.position = "none",
        legend.text=element_text(size=9),
        legend.title=element_text(size=9)) +
  scale_color_manual(values=c("grey47", "black")) +
  scale_y_continuous(name="", limits=c(-1.2,1))  



# 2012
p_natid_m12<- ggplot(natid_m12, aes(color=NatID, group=NatID)) +  ggtitle("2012") +
  geom_pointrange(aes(x = Variable, y = Coefficient,  group=NatID, 
                      ymin = Coefficient - SE*interval, ymax = Coefficient + SE*interval),
                  position = position_dodge(width = 1/1.6), shape=rep(c(15,17), each=5)) +
  coord_flip() + theme_classic() + labs(y="", x="") +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_vline(xintercept = 2.5, color = gray(1/2), lty = 2) +
  guides(color=guide_legend(title="(ref: neutral)")) +
  theme(plot.title = element_text(size = 10),
        legend.position = "none",
        legend.text=element_text(size=9),
        legend.title=element_text(size=9)) +
  scale_color_manual(values=c("grey47", "black")) +
  scale_y_continuous(name="", limits=c(-1.2,1))  



# 2009
p_natid_m09<- ggplot(natid_m09, aes(color=NatID, group=NatID)) +  ggtitle("2009") +
  geom_pointrange(aes(x = Variable, y = Coefficient,  group=NatID, 
                      ymin = Coefficient - SE*interval, ymax = Coefficient + SE*interval),
                  position = position_dodge(width = 1/1.6), shape=rep(c(15,17), each=5)) +
  coord_flip() + theme_classic() + labs(y="", x="") +
  geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
  geom_vline(xintercept = 2.5, color = gray(1/2), lty = 2) +
  guides(color=guide_legend(title="(ref: neutral)")) +
  theme(plot.title = element_text(size = 10),
        legend.position = "none",
        legend.text=element_text(size=9),
        legend.title=element_text(size=9)) +
  scale_color_manual(values=c("grey47", "black")) +
  scale_y_continuous(name="", limits=c(-1.2,1)) 


#---------------------------
# Combine Plots and export: 
#---------------------------
pdf(file="figureA1.pdf", width=6, height=5)  
grid.arrange(p_natid_m20, p_natid_m16,p_natid_m12, p_natid_m09,  nrow = 2)
dev.off()



#===================================================
# Figure A2: Model Estimates on Ideological Groups
#===================================================

# 2020

mat_ideoself_m20 <- matrix(NA, 20, ncol = 4)
colnames(mat_ideoself_m20) <- c("Variable","Coefficient", "SE", "IdeoSelf")
mat_ideoself_m20[,1] <- c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
mat_ideoself_m20[,2] <- t(coef.m20[6:9,1:5])
mat_ideoself_m20[,3] <- t(se_m20[6:9,1:5]) 
mat_ideoself_m20[,4] <- rep(c("Conservative", "Socialist", "Feminist.Ecologist", "Liberal"), each=5) 
mat_ideoself_m20

ideoself_m20<-as.data.frame(mat_ideoself_m20)
ideoself_m20$Coefficient <-as.numeric(ideoself_m20$Coefficient)
ideoself_m20$SE <-as.numeric(ideoself_m20$SE)
ideoself_m20
ideoself_m20$Variable <- factor(ideoself_m20$Variable, c("PNV", "EH-Bildu", "PP", "PSE",  "Podemos"))

# 2016

mat_ideoself_m16 <- matrix(NA, 20, ncol = 4)
colnames(mat_ideoself_m16) <- c("Variable","Coefficient", "SE", "IdeoSelf")
mat_ideoself_m16[,1] <- c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
mat_ideoself_m16[,2] <- t(coef.m16[6:9,1:5])
mat_ideoself_m16[,3] <- t(se_m16[6:9,1:5]) 
mat_ideoself_m16[,4] <- rep(c("Conservative", "Socialist", "Feminist.Ecologist", "Liberal"), each=5) 
mat_ideoself_m16

ideoself_m16<-as.data.frame(mat_ideoself_m16)
ideoself_m16$Coefficient <-as.numeric(ideoself_m16$Coefficient)
ideoself_m16$SE <-as.numeric(ideoself_m16$SE)
ideoself_m16
ideoself_m16$Variable <- factor(ideoself_m16$Variable, c("PNV", "EH-Bildu", "PP", "PSE",  "Podemos"))


# 2012

mat_ideoself_m12 <- matrix(NA, 20, ncol = 4)
colnames(mat_ideoself_m12) <- c("Variable","Coefficient", "SE", "IdeoSelf")
mat_ideoself_m12[,1] <- c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
mat_ideoself_m12[,2] <- t(coef.m12[6:9,1:5])
mat_ideoself_m12[,3] <- t(se_m12[6:9,1:5]) 
mat_ideoself_m12[,4] <- rep(c("Conservative", "Socialist", "Feminist.Ecologist", "Liberal"), each=5) 
mat_ideoself_m12

ideoself_m12<-as.data.frame(mat_ideoself_m12)
ideoself_m12$Coefficient <-as.numeric(ideoself_m12$Coefficient)
ideoself_m12$SE <-as.numeric(ideoself_m12$SE)
ideoself_m12
ideoself_m12$Variable <- factor(ideoself_m12$Variable, c("PNV", "EH-Bildu", "PP", "PSE",  "Podemos"))


# 2009

mat_ideoself_m09 <- matrix(NA, 20, ncol = 4)
colnames(mat_ideoself_m09) <- c("Variable","Coefficient", "SE", "IdeoSelf")
mat_ideoself_m09[,1] <- c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
mat_ideoself_m09[,2] <- t(coef.m09[6:9,1:5])
mat_ideoself_m09[,3] <- t(se_m09[6:9,1:5]) 
mat_ideoself_m09[,4] <- rep(c("Conservative", "Socialist", "Feminist.Ecologist", "Liberal"), each=5) 
mat_ideoself_m09


ideoself_m09<-as.data.frame(mat_ideoself_m09)
ideoself_m09$Coefficient <-as.numeric(ideoself_m09$Coefficient)
ideoself_m09$SE <-as.numeric(ideoself_m09$SE)
ideoself_m09
ideoself_m09$Variable <- factor(ideoself_m09$Variable, c("PNV", "EH-Bildu", "PP", "PSE",  "Podemos"))


#---------------------
# create plots
#---------------------
# 2020
p_ideoself_m20<- ggplot(ideoself_m20, aes(color=IdeoSelf, group=IdeoSelf)) + ggtitle("2020") + 
   geom_pointrange(aes(x = Variable, y = Coefficient,  group=IdeoSelf, 
                       ymin = Coefficient - SE*interval, ymax = Coefficient + SE*interval),
                   position = position_dodge(width = 1/1.3), shape=rep(c(15,1,17,5), each=5)) +
   coord_flip() + theme_classic() + labs(y="", x="") +
   geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
   geom_vline(xintercept = 2.5, color = gray(1/2), lty = 2) +
   #guides(color=guide_legend(title="", override.aes=list(shape =c(15,1,17,5)))) +
   theme(plot.title = element_text(size = 10),
         legend.position = "none",
         legend.text=element_text(size=8),
         legend.title=element_text(size=8)) +
   scale_color_manual(values=c("gray35", "gray35","gray35", "gray35" )) +
   scale_y_continuous(name="", limits=c(-1.2,1.6), breaks=c(-1,-0.5,0,0.5,1,1.5)) 


# 2016
p_ideoself_m16<- ggplot(ideoself_m16, aes(color=IdeoSelf, group=IdeoSelf)) +  ggtitle("2016") +
   geom_pointrange(aes(x = Variable, y = Coefficient,  group=IdeoSelf, 
                       ymin = Coefficient - SE*interval, ymax = Coefficient + SE*interval),
                   position = position_dodge(width = 1/1.3), shape=rep(c(15,1,17,5), each=5)) +
   coord_flip() + theme_classic() + labs(y="", x="") +
   geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
   geom_vline(xintercept = 2.5, color = gray(1/2), lty = 2) +
   guides(color=guide_legend(title="", override.aes=list(shape =c(15,1,17,5)))) +  
   theme(plot.title = element_text(size = 10),
         legend.position = "none",
         legend.text=element_text(size=8),
         legend.title=element_text(size=8)) +
   scale_color_manual(values=c("gray35", "gray35","gray35", "gray35" )) +
   scale_y_continuous(name="", limits=c(-1.2,1.6), breaks=c(-1,-0.5,0,0.5,1,1.5))  


# 2012
p_ideoself_m12<- ggplot(ideoself_m12, aes(color=IdeoSelf, group=IdeoSelf)) +  ggtitle("2012") +
   geom_pointrange(aes(x = Variable, y = Coefficient,  group=IdeoSelf, 
                       ymin = Coefficient - SE*interval, ymax = Coefficient + SE*interval),
                   position = position_dodge(width = 1/1.3), shape=rep(c(15,1,17,5), each=5)) +
   coord_flip() + theme_classic() + labs(y="", x="") +
   geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
   geom_vline(xintercept = 2.5, color = gray(1/2), lty = 2) +
   guides(color=guide_legend(title="(ref: Nationalist)", override.aes=list(shape =c(15,1,17,5)))) +
   theme(plot.title = element_text(size = 10),
         legend.position = "none",
         legend.text=element_text(size=8),
         legend.title=element_text(size=8)) +
   scale_color_manual(values=c("gray35", "gray35","gray35", "gray35" )) +
   scale_y_continuous(name="", limits=c(-1.2,1.6), breaks=c(-1,-0.5,0,0.5,1,1.5)) 



# 2009
p_ideoself_m09<- ggplot(ideoself_m09, aes(color=IdeoSelf, group=IdeoSelf)) +  ggtitle("2009") +
   geom_pointrange(aes(x = Variable, y = Coefficient,  group=IdeoSelf, 
                       ymin = Coefficient - SE*interval, ymax = Coefficient + SE*interval),
                   position = position_dodge(width = 1/1.3), shape=rep(c(15,1,17,5), each=5)) +
   coord_flip() + theme_classic() + labs(y="", x="") +
   geom_hline(yintercept = 0, color = gray(1/2), lty = 2) +
   geom_vline(xintercept = 2.5, color = gray(1/2), lty = 2) +
   guides(color=guide_legend(title="(ref: Nationalist)", override.aes=list(shape =c(15,1,17,5)))) +
   theme(plot.title = element_text(size = 10),
         legend.position = "none",
         legend.text=element_text(size=8),
         legend.title=element_text(size=8 )) +
   scale_color_manual(values=c("gray35", "gray35","gray35", "gray35" )) +
   scale_y_continuous(name="", limits=c(-1.2,1.6), breaks=c(-1,-0.5,0,0.5,1,1.5))  


#---------------------------
# Combine Plots and export: 
#---------------------------
pdf(file="figureA2.pdf", width=6, height=7)  
grid.arrange(p_ideoself_m20, p_ideoself_m16, p_ideoself_m12, p_ideoself_m09,  nrow = 2)
dev.off()


#=====================================================
# Figure A3: Model Estimates on Catholic Denomination
#=====================================================

# 2020

mat_cath_m20 <- matrix(NA, 5, ncol = 3)
colnames(mat_cath_m20) <- c("Variable","Coefficient", "SE")
mat_cath_m20[,1] <- c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
mat_cath_m20[,2] <- t(coef.m20[15:15,1:5])
mat_cath_m20[,3] <- t(se_m20[15:15,1:5]) 
mat_cath_m20

cath_m20<-as.data.frame(mat_cath_m20)
cath_m20$Coefficient <-as.numeric(cath_m20$Coefficient)
cath_m20$SE <-as.numeric(cath_m20$SE)
cath_m20
cath_m20$Variable <- factor(cath_m20$Variable, c("PNV", "EH-Bildu", "PP", "PSE",  "Podemos"))


# 2016

mat_cath_m16 <- matrix(NA, 5, ncol = 3)
colnames(mat_cath_m16) <- c("Variable","Coefficient", "SE")
mat_cath_m16[,1] <- c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
mat_cath_m16[,2] <- t(coef.m16[15:15,1:5])
mat_cath_m16[,3] <- t(se_m16[15:15,1:5]) 
mat_cath_m16

cath_m16<-as.data.frame(mat_cath_m16)
cath_m16$Coefficient <-as.numeric(cath_m16$Coefficient)
cath_m16$SE <-as.numeric(cath_m16$SE)
cath_m16
cath_m16$Variable <- factor(cath_m16$Variable, c("PNV", "EH-Bildu", "PP", "PSE",  "Podemos"))


# 2012

mat_cath_m12 <- matrix(NA, 5, ncol = 3)
colnames(mat_cath_m12) <- c("Variable","Coefficient", "SE")
mat_cath_m12[,1] <- c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
mat_cath_m12[,2] <- t(coef.m12[15:15,1:5])
mat_cath_m12[,3] <- t(se_m12[15:15,1:5]) 
mat_cath_m12

cath_m12<-as.data.frame(mat_cath_m12)
cath_m12$Coefficient <-as.numeric(cath_m12$Coefficient)
cath_m12$SE <-as.numeric(cath_m12$SE)
cath_m12
cath_m12$Variable <- factor(cath_m12$Variable, c("PNV", "EH-Bildu", "PP", "PSE",  "Podemos"))



#---------------------
# create plots
#---------------------

# 2020
p_cath_m20 <- ggplot(cath_m20) +  
   ggtitle("2020") + 
   geom_pointrange(aes(x = Variable, y = Coefficient,  ymin = Coefficient - SE*interval, ymax = Coefficient + SE*interval),
                   lwd = 1/2, position = position_dodge(width = 1/2)) +    
   geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + coord_flip() + theme_classic() + labs(x="", y="") +
   geom_vline(xintercept = 2.5, color = gray(1/2), lty = 2) +
   theme(legend.position = "top") + theme(plot.title = element_text(size = 10)) +
   scale_y_continuous(name="", limits=c(-1,0.6), breaks=c(-1,-0.5,0,0.5))


# 2016
p_cath_m16 <- ggplot(cath_m16) +  
   ggtitle("2016") + 
   geom_pointrange(aes(x = Variable, y = Coefficient,  ymin = Coefficient - SE*interval, ymax = Coefficient + SE*interval),
                   lwd = 1/2, position = position_dodge(width = 1/2)) +    
   geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + coord_flip() + theme_classic() + labs(x="", y="") +
   geom_vline(xintercept = 2.5, color = gray(1/2), lty = 2) +
   theme(legend.position = "top") + theme(plot.title = element_text(size = 10)) +
   scale_y_continuous(name="", limits=c(-1,0.6), breaks=c(-1,-0.5,0,0.5))

# 2012
p_cath_m12 <- ggplot(cath_m12) +  
   ggtitle("2012") + 
   geom_pointrange(aes(x = Variable, y = Coefficient,  ymin = Coefficient - SE*interval, ymax = Coefficient + SE*interval),
                   lwd = 1/2, position = position_dodge(width = 1/2)) +    
   geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + coord_flip() + theme_classic() + labs(x="", y="") +
   geom_vline(xintercept = 2.5, color = gray(1/2), lty = 2) +
   theme(legend.position = "top") + theme(plot.title = element_text(size = 10)) +
   scale_y_continuous(name="", limits=c(-1,0.6), breaks=c(-1,-0.5,0,0.5))


#---------------------------
# Combine Plots and export: 
#---------------------------
pdf(file="figureA3.pdf", width=9, height=2.5)  
grid.arrange(p_cath_m20, p_cath_m16, p_cath_m12, nrow = 1)
dev.off()


#====================================================================
# Figure A4: Model Estimates on Basque Language Usage (additive index)
#=====================================================================

# 2020

mat_eusk_m20 <- matrix(NA, 5, ncol = 3)
colnames(mat_eusk_m20) <- c("Variable","Coefficient", "SE")
mat_eusk_m20[,1] <- c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
mat_eusk_m20[,2] <- t(coef.m20[4:4,1:5])
mat_eusk_m20[,3] <- t(se_m20[4:4,1:5]) 
mat_eusk_m20

eusk_m20<-as.data.frame(mat_eusk_m20)
eusk_m20$Coefficient <-as.numeric(eusk_m20$Coefficient)
eusk_m20$SE <-as.numeric(eusk_m20$SE)
eusk_m20
eusk_m20$Variable <- factor(eusk_m20$Variable, c("PNV", "EH-Bildu", "PP", "PSE",  "Podemos"))


# 2016

mat_eusk_m16 <- matrix(NA, 5, ncol = 3)
colnames(mat_eusk_m16) <- c("Variable","Coefficient", "SE")
mat_eusk_m16[,1] <- c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
mat_eusk_m16[,2] <- t(coef.m16[4:4,1:5])
mat_eusk_m16[,3] <- t(se_m16[4:4,1:5]) 
mat_eusk_m16

# save as data frame
eusk_m16<-as.data.frame(mat_eusk_m16)
eusk_m16$Coefficient <-as.numeric(eusk_m16$Coefficient)
eusk_m16$SE <-as.numeric(eusk_m16$SE)
eusk_m16
eusk_m16$Variable <- factor(eusk_m16$Variable, c("PNV", "EH-Bildu", "PP", "PSE",  "Podemos"))


# 2012

mat_eusk_m12 <- matrix(NA, 5, ncol = 3)
colnames(mat_eusk_m12) <- c("Variable","Coefficient", "SE")
mat_eusk_m12[,1] <- c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
mat_eusk_m12[,2] <- t(coef.m12[4:4,1:5])
mat_eusk_m12[,3] <- t(se_m12[4:4,1:5]) 
mat_eusk_m12

eusk_m12<-as.data.frame(mat_eusk_m12)
eusk_m12$Coefficient <-as.numeric(eusk_m12$Coefficient)
eusk_m12$SE <-as.numeric(eusk_m12$SE)
eusk_m12
eusk_m12$Variable <- factor(eusk_m12$Variable, c("PNV", "EH-Bildu", "PP", "PSE",  "Podemos"))

# 2009

mat_eusk_m09 <- matrix(NA, 5, ncol = 3)
colnames(mat_eusk_m09) <- c("Variable","Coefficient", "SE")
mat_eusk_m09[,1] <- c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
mat_eusk_m09[,2] <- t(coef.m09[4:4,1:5])
mat_eusk_m09[,3] <- t(se_m09[4:4,1:5]) 
mat_eusk_m09

eusk_m09<-as.data.frame(mat_eusk_m09)
eusk_m09$Coefficient <-as.numeric(eusk_m09$Coefficient)
eusk_m09$SE <-as.numeric(eusk_m09$SE)
eusk_m09
eusk_m09$Variable <- factor(eusk_m09$Variable, c("PNV", "EH-Bildu", "PP", "PSE",  "Podemos"))


#---------------------
# create plots
#---------------------

# 2020
p_eusk_m20 <- ggplot(eusk_m20) +  
   ggtitle("2020") + 
   geom_pointrange(aes(x = Variable, y = Coefficient,  ymin = Coefficient - SE*interval, ymax = Coefficient + SE*interval),
                   lwd = 1/2, position = position_dodge(width = 1/2), color="grey47", shape=17) +    
   geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + coord_flip() + theme_classic() + labs( x="", y="") +
   geom_vline(xintercept = 2.5, color = gray(1/2), lty = 2) +
   theme(legend.position = "top") + theme(plot.title = element_text(size = 10)) +
   scale_y_continuous(name="", limits=c(-1,1.2), breaks=c(-1,-0.5,0,0.5, 1))


# 2016
p_eusk_m16 <- ggplot(eusk_m16) +  
   ggtitle("2016") + 
   geom_pointrange(aes(x = Variable, y = Coefficient,  ymin = Coefficient - SE*interval, ymax = Coefficient + SE*interval),
                   lwd = 1/2, position = position_dodge(width = 1/2), color="grey47", shape=17) +    
   geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + coord_flip() + theme_classic() + labs(x="", y="") +
   geom_vline(xintercept = 2.5, color = gray(1/2), lty = 2) +
   theme(legend.position = "top") + theme(plot.title = element_text(size = 10)) +
   scale_y_continuous(name="", limits=c(-1,1.2), breaks=c(-1,-0.5,0,0.5, 1))

# 2012
p_eusk_m12 <- ggplot(eusk_m12) +  
   ggtitle("2012") + 
   geom_pointrange(aes(x = Variable, y = Coefficient,  ymin = Coefficient - SE*interval, ymax = Coefficient + SE*interval),
                   lwd = 1/2, position = position_dodge(width = 1/2), color="grey47", shape=17) +    
   geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + coord_flip() + theme_classic() + labs(x="", y="") +
   geom_vline(xintercept = 2.5, color = gray(1/2), lty = 2) +
   theme(legend.position = "top") + theme(plot.title = element_text(size = 10)) +
   scale_y_continuous(name="", limits=c(-1,1.2), breaks=c(-1,-0.5,0,0.5, 1))

# 2009
p_eusk_m09 <- ggplot(eusk_m09) +  
   ggtitle("2009") + 
   geom_pointrange(aes(x = Variable, y = Coefficient,  ymin = Coefficient - SE*interval, ymax = Coefficient + SE*interval),
                   lwd = 1/2, position = position_dodge(width = 1/2), color="grey47", shape=17) +    
   geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + coord_flip() + theme_classic() + labs(x="", y="") +
   geom_vline(xintercept = 2.5, color = gray(1/2), lty = 2) +
   theme(legend.position = "top") + theme(plot.title = element_text(size = 10)) +
   scale_y_continuous(name="", limits=c(-1,1.2), breaks=c(-1,-0.5,0,0.5, 1))


#---------------------------
# Combine Plots and export: 
#---------------------------
pdf(file="figureA4.pdf", width=6, height=5)  
grid.arrange(p_eusk_m20, p_eusk_m16, p_eusk_m12, p_eusk_m09, nrow = 2)
dev.off()


#============================================================
# Figure A5: Model Estimates on Left-Right Self-Placements
#============================================================

# 2020

mat_lrego_m20 <- matrix(NA, 5, ncol = 3)
colnames(mat_lrego_m20) <- c("Variable","Coefficient", "SE")
mat_lrego_m20[,1] <- c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
mat_lrego_m20[,2] <- t(coef.m20[5:5,1:5])
mat_lrego_m20[,3] <- t(se_m20[5:5,1:5]) 
mat_lrego_m20

lrego_m20<-as.data.frame(mat_lrego_m20)
lrego_m20$Coefficient <-as.numeric(lrego_m20$Coefficient)
lrego_m20$SE <-as.numeric(lrego_m20$SE)
lrego_m20
lrego_m20$Variable <- factor(lrego_m20$Variable, c("PNV", "EH-Bildu", "PP", "PSE",  "Podemos"))


# 2016

mat_lrego_m16 <- matrix(NA, 5, ncol = 3)
colnames(mat_lrego_m16) <- c("Variable","Coefficient", "SE")
mat_lrego_m16[,1] <- c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
mat_lrego_m16[,2] <- t(coef.m16[5:5,1:5])
mat_lrego_m16[,3] <- t(se_m16[5:5,1:5]) 
mat_lrego_m16

lrego_m16<-as.data.frame(mat_lrego_m16)
lrego_m16$Coefficient <-as.numeric(lrego_m16$Coefficient)
lrego_m16$SE <-as.numeric(lrego_m16$SE)
lrego_m16
lrego_m16$Variable <- factor(lrego_m16$Variable, c("PNV", "EH-Bildu", "PP", "PSE",  "Podemos"))


# 2012

mat_lrego_m12 <- matrix(NA, 5, ncol = 3)
colnames(mat_lrego_m12) <- c("Variable","Coefficient", "SE")
mat_lrego_m12[,1] <- c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
mat_lrego_m12[,2] <- t(coef.m12[5:5,1:5])
mat_lrego_m12[,3] <- t(se_m12[5:5,1:5]) 
mat_lrego_m12

lrego_m12<-as.data.frame(mat_lrego_m12)
lrego_m12$Coefficient <-as.numeric(lrego_m12$Coefficient)
lrego_m12$SE <-as.numeric(lrego_m12$SE)
lrego_m12
lrego_m12$Variable <- factor(lrego_m12$Variable, c("PNV", "EH-Bildu", "PP", "PSE",  "Podemos"))

# 2009

mat_lrego_m09 <- matrix(NA, 5, ncol = 3)
colnames(mat_lrego_m09) <- c("Variable","Coefficient", "SE")
mat_lrego_m09[,1] <- c("Podemos", "PSE", "PP", "EH-Bildu", "PNV")
mat_lrego_m09[,2] <- t(coef.m09[5:5,1:5])
mat_lrego_m09[,3] <- t(se_m09[5:5,1:5]) 
mat_lrego_m09

lrego_m09<-as.data.frame(mat_lrego_m09)
lrego_m09$Coefficient <-as.numeric(lrego_m09$Coefficient)
lrego_m09$SE <-as.numeric(lrego_m09$SE)
lrego_m09
lrego_m09$Variable <- factor(lrego_m09$Variable, c("PNV", "EH-Bildu", "PP", "PSE",  "Podemos"))

#---------------------
# create plots
#---------------------

# 2020
p_lrego_m20 <- ggplot(lrego_m20) +  
  ggtitle("2020") + 
  geom_pointrange(aes(x = Variable, y = Coefficient,  ymin = Coefficient - SE*interval, ymax = Coefficient + SE*interval),
                  lwd = 1/2, position = position_dodge(width = 1/2)) +    
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + coord_flip() + theme_classic() + labs(x="", y="") +
  geom_vline(xintercept = 2.5, color = gray(1/2), lty = 2) +
  theme(legend.position = "top") + theme(plot.title = element_text(size = 10)) +
  scale_y_continuous(name="", limits=c(-0.4,0.4), breaks=c(-0.4,-0.2,0,0.2,0.4))


# 2016
p_lrego_m16 <- ggplot(lrego_m16) +  
  ggtitle("2016") + 
  geom_pointrange(aes(x = Variable, y = Coefficient,  ymin = Coefficient - SE*interval, ymax = Coefficient + SE*interval),
                  lwd = 1/2, position = position_dodge(width = 1/2)) +    
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + coord_flip() + theme_classic() + labs(x="", y="") +
  geom_vline(xintercept = 2.5, color = gray(1/2), lty = 2) +
  theme(legend.position = "top") + theme(plot.title = element_text(size = 10)) +
  scale_y_continuous(name="", limits=c(-0.4,0.4), breaks=c(-0.4,-0.2,0,0.2,0.4))

# 2012
p_lrego_m12 <- ggplot(lrego_m12) +  
  ggtitle("2012") + 
  geom_pointrange(aes(x = Variable, y = Coefficient,  ymin = Coefficient - SE*interval, ymax = Coefficient + SE*interval),
                  lwd = 1/2, position = position_dodge(width = 1/2)) +    
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + coord_flip() + theme_classic() + labs(x="", y="") +
  geom_vline(xintercept = 2.5, color = gray(1/2), lty = 2) +
  theme(legend.position = "top") + theme(plot.title = element_text(size = 10)) +
  scale_y_continuous(name="", limits=c(-0.4,0.4), breaks=c(-0.4,-0.2,0,0.2,0.4))

# 2009
p_lrego_m09 <- ggplot(lrego_m09) +  
  ggtitle("2009") + 
  geom_pointrange(aes(x = Variable, y = Coefficient,  ymin = Coefficient - SE*interval, ymax = Coefficient + SE*interval),
                  lwd = 1/2, position = position_dodge(width = 1/2)) +    
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + coord_flip() + theme_classic() + labs(x="", y="") +
  geom_vline(xintercept = 2.5, color = gray(1/2), lty = 2) +
  theme(legend.position = "top") + theme(plot.title = element_text(size = 10)) +
  scale_y_continuous(name="", limits=c(-0.4,0.4), breaks=c(-0.4,-0.2,0,0.2,0.4))

#---------------------------
# Combine Plots and export: 
#---------------------------
pdf(file="figureA5.pdf", width=6, height=5)  
grid.arrange(p_lrego_m20, p_lrego_m16, p_lrego_m12, p_lrego_m09,  nrow = 2)
dev.off()

